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We show that the Riccati form of the one-dimensional Schrodinger equation can be reformulated 
in terms of two linear equations depending on an arbitrary function G. When G and the potential (as 
for anharmonic oscillators) are polynomials the solutions of these two equations are entire functions 
(L and K) and the zeroes of K are identical to those of the wave function. Requiring such a zero at 
a large but finite value of the argument yields the low energy eigenstates with exponentially small 
errors. Approximate formulas for these errors are provided. We explain how to chose G in order 
improve dramatically the numerical treatment. The method yields many significant digits with 
modest computer means. We discuss the extension of this method in the case of several variables. 

PACS: 03.65.Ge, 03.65.-w, 02.30.Em, 02.30.Mv, 10.10.St, 33.20.-t 



I. INTRODUCTION 



Quantum anharmonic oscillators appear in a wide vari- 
ety of problems in molecular, nuclear or condensed mat- 
ter physics. Typically, anharmonic terms appear in ex- 
pansions about a minimum of a potential, when ones tries 
to incorporate the non-linear features of the forces re- 
sponsible for this equilibrium. The most celebrated ex- 
ample is the quartic anharmonic oscillator [1] where a 
Ax 4 term is added to the usual harmonic Hamiltonian. 
Introducing bilinear couplings among a set of such oscil- 
lators leads to a rich spectrum, for instance, multiphonon 
bound states in one-dimensional lattice models [2] . More 
generally, one can think about the A</> 4 (or higher powers 
of <j)) field theories in various dimensions as systems of 
coupled anharmonic oscillators. 

Anharmonic terms can be treated perturbatively and 
the perturbative series can be represented by Feynman 
diagrams. Unfortunately, the coefficients of the series 
[1,3] have a factorial growth and the numerical values ob- 
tained from the truncated series have an accuracy which 
is subject to limitations. At fixed coupling, there is an 
order at which an optimal accuracy is reached. At fixed 
order, there is a value of the coupling beyond which the 
numerical values are meaningless even as an order of mag- 
nitude. In the case of the single-well quartic potential, 
Pade approximants can be used for the series or its Borel 
transform. Rigorous proofs of convergence can be es- 
tablished in particular cases [4]. Unfortunately, such a 
method does not apply to the case of the double-well po- 
tential [5] where instanton effects [6,7] need to be taken 
into account. It should also be noted that even when 
Pade approximants converge, the convergence rate may 
be slow. Strong coupling expansions [8] or variational in- 
terpolations [9] sometimes provide more accurate results. 

The above discussion shows that finding an expansion 
which can be used indiscriminately for most quantum 
mechanical problems with polynomial potentials remains 
a challenging problem. Alternatively, one can use nu- 
merical methods. Variational methods are often used to 
obtain upper and lower bounds on energy levels [10,11]. 



These methods are based on rigorous inequalities and are 
considered superior to methods based on numerical inte- 
gration [11]. However, the difference between the bounds 
widens rapidly with the anharmonic coupling and the en- 
ergy level. Methods based on series expansions in the 
position variable [12-15] appear to produce more signif- 
icant digits more easily. However, our understanding of 
the convergence and numerical stability of these methods 
seems to be limited to empirical observations. The meth- 
ods based on series expansions fall into two categories: 
methods based on the evaluations of determinants [12,14] 
and methods based on boundary conditions at large but 
finite values of the position [13,15] . The main goal of this 
article is to provide a systematic discussion of the errors 
associated with this second category of methods and to 
show how to make these errors arbitrarily small in the 
most efficient way. With the exception of Section IX, we 
only consider one-dimensional problems. We discuss two 
types of errors. First, the numerical errors made in calcu- 
lating the energy which makes the wave function vanish 
at some large value of the position x max . Second, the 
intrinsic error due to the fmiteness of x max . 

The basic elements the numerical method used here- 
after were sketched in Ref. [15] and applied to the quartic 
anharmonic oscillator. We wrote the logarithmic deriva- 
tive of the wave function which appears in the Riccati 
equation as L/K and showed that these functions were 
entire. The values of the first ten eigenvalues with 30 
significant digits provided for a particular coupling have 
been used to test new theoretical methods [16]. Two is- 
sues were left open in this formulation: first, the basic 
equations had an interesting invariance which was not 
undestood but could be used to improve the numerical 
efficiency; second, the use of the method for parity non- 
invariant potentials appeared to be unduly complicated 
[17]. 

In Section II, we present a new formulation where these 
two issues are settled. The basic equations presented de- 
pend on an arbitrary function denoted G(x). This free- 
dom can be interpreted as a local gauge invariance as- 
sociated with the fact that only L/K is physical. The 
wave function is invariant under these local transforma- 
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tions. In section III, we show how to construct power 
series for L and K. The complications in the case of par- 
ity non-invariant potentials (such as asymmetric double- 
wells) are minimal. When the potential and the gauge 
function are polynomials, these series define entire func- 
tion. In other words, it is always possible to construct 
arbitrarily accurate solutions of the Schrodinger equation 
for arbitrary E within a given range of the position vari- 
able, by calculating enough terms in the expansions of 
L and K. This allows us to reproduce the asymptotic 
behavior of the wave function and determine the energy 
eigenvalues. In section IV, we use the global properties 
of the flows of the Riccati equation to recall of some ba- 
sic results related to the WKB approximation and the 
Sturm-Liouville theorem. We explain how bifurcations 
in the asymptotic behavior of the functions K and L can 
be exploited to determine the eigenvalues. 

It should be noted that the importance of reproducing 
the proper asymptotic behavior has been emphasized in 
variational approaches [18]. It should also be noted that 
Pade approximants have been used in conjunction with 
the Riccati equation in Ref. [14], where the quantization 
condition used was that the approximants give one addi- 
tional coefficient in the Taylor expansion. This procedure 
depends only on the coefficients of the expansions used 
and there is no reference to any particular value of x (as 
our x max ). Consequently, there is no obvious connection 
between the two approaches. 

In the next two sections, we show how to turn the 
gauge invariance to our advantage. In Section V, the 
quantitative aspects of the bifurcation are discussed with 
an exponential parametrization similar to the one used to 
determine Lyapounov exponents in the study of chaotic 
dynamical system. The exponents are G-dependent. We 
provide an approximate way to determine the exponents 
and the energy resolution. We explain how our freedom 
in chosing G can be used to make the bifurcation more 
violent and improve the energy resolution. However, the 
choice of G also affects the convergence of L and K and 
consequently the numerical accuracy of the solution of 
the Schrodinger equation. In Section VI, we show in a 
particular example that for an expansion of L and K at 
a given order, a judicious choice of gauge can improve 
tremendously the numerical accuracy of an energy level. 
We discuss the two principles which allow to make opti- 
mal choices of G and provide practical methods to deter- 
mine approximately this optimal choice for the general 
case. We use these methods to explain some empirical 
results found in [13]. 

In Section VII, we discuss the the error SE on the en- 
ergy levels due to the finiteness of x max . We propose two 
approximate formulas valid, respectively, for intermedi- 
ate and large values of x max and compatible in overlap- 
ping ranges. Note that one can reinterpret the condi- 
tion that the wavefunction vanishes at x max as coming 
from a slightly different problem where the potential be- 
comes infinite at x max . In the path- integral formulation 
(which can be extended immediately to field theory prob- 



lems), the fact that the potential becomes infinite at x max 
means that paths with values of x larger than x max are 
not taken into account. It has been argued [19,20] that 
these configurations are responsible for the asymptotic 
behavior of the regular perturbative series. In Ref. [20], 
we showed that the perturbative series of several mod- 
ified problem were convergent. The error formula sets 
the accuracy limitations of this approach. Some of the 
methods used in this section could be used for quantum 
field theory problems. 

The anharmonic oscillator can be considered as a field 
theory with one time and zero space dimensions. It can 
be used to test approximate methods such as perturba- 
tive expansions or semi-classical procedures. An illustra- 
tive example is given in Ref. [23] where multi-instanton 
effects were considered and where the splitting of the 
two lowest levels of a double- well problem were estimated 
with more than hundred digits. In Section VIII, we show 
that our method can be used to reproduce all these digits. 
Finally, we discuss the generalization of the method to 
problems with several variables in Section IX. For these 
problems, our ability to reduce the degree of expansion 
by using optimal gauge functions may be crucial. 



II. BASIC EQUATIONS AND THEIR 
GAUGE-INVARIANCE 

We consider a one-dimensional, time-independent 
Schrodinger equation H^> = E^, for an Hamiltonian 

2 21 

H=^ + Y [ V 3 x>. (1) 

As is well-known, one can reexpress the wave function in 
terms of its logarithmic derivative 

*(z) oc e^> 0fe) ' (2) 

and obtain the Riccati form of the equation: 

h<t>' = 4> 2 + 2m(E - V) . (3) 

It is assumed that m > and that the leading power of 
V is even with a positive coefficient (V21 > 0). 

Writing <f> = L/K, we obtain a solution of Eq. (3) 
provided that we solve the system of equations: 

TiL' + 2m(V - E)K + GL = (4) 
fiK' + L + GK = (5) 

where G{x) is an unspecified function. This can be seen 
by multiplying (4) by K, (5) by L and eliminating GKL 
by taking the difference. One then obtains the Riccati 
equation (3) multiplied by K 2 . Near a zero of K, one 
can check that Eqs. (4-5) remain valid, namely they im- 
pose that cj> has a simple pole with residue — K. This 
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allows the wave function to become zero and change sign 
as the contour goes around the pole on either side. 
Eqs. (4-5) are invariant under the local transformation 



K n+ i — 



-1 



L(x) -» Q{x)L{x) 
K(x) -» Q(x)K(x) 
G{x)^G{x)-hQ'{x)/Q{x) 



(6) 



where Q(x) is an arbitrary function. It is clear that 
this transformation leaves 4> and the wave function un- 
changed. If we choose G = and eliminate L using Eq. 
(5), we recover the Schrodingcr equation for K. Starting 
from this gauge and making an arbitrary transformation, 
we find that in general 



K{x) oc *(x)e h Jx ° 



i P dyG(y) 



(7) 



This shows that when G is polynomial, K is simply ^ 
multiplied by an entire function with no zeroes [21]. This 
means that the zeroes of K and \t are identical. In other 
words, there are no spurious zeroes when G is polyno- 
mial. 

By taking the derivative of Eqs. (4) and (5) and choos- 
ing G{x) appropriately, one can obtain the basic Equa- 
tions used in [15]. The explicit form of G{x) is reached by 
comparing the two sets of equations and integrating one 
of the differences. The two possibilities are compatible. 
The resulting integral expression can be worked out eas- 
ily by the interested reader. The only important point is 
that the G found that way is in general not polynomial, 
justifying the spurious zeroes found with the original for- 
mulation. 



III. SOLUTIONS IN TERMS OF ENTIRE 
FUNCTIONS 

The function G can be chosen at our convenience. For 
instance, we could impose the condition K = 1 by taking 
G = —L and recover the Riccati equation for L. How- 
ever, the main advantage of Eqs. (4-5) is that they are 
linear first order differential equations with variables co- 
efficients. It is well-known [22] that if we consider these 
equations for complex x, the solutions inherit the do- 
main of analyticity of the coefficients (provided that this 
domain is simply connected). If the coefficients are en- 
tire functions, there exists a unique entire solution cor- 
responding to a particular set of initial values. In the 
following, we restrict ourselves to the case where V and 
G are polynomials. 

One can construct the unique solution corresponding 
to a particular choice of initial values L(0) and K(0) by 
series expansions. Using K(x) — Y^=o K n x n and similar 
notations for the other functions, one obtains the simple 
recursion 



L n +i — 



-1 



h(n+l) 



-( (2mV l K p + L l G p )-2mEK n ) 



h(n+l) 



(L n + K i°p) 

l-\-p—n 



(8) 



Given Lq and Kq, these equations allow to determine all 
the other coefficients. For potentials which are parity in- 
variant, and if G is an odd function, L and K can be 
assigned definite and opposite parities. In this case, we 
can impose the initial conditions Kq — 1 and Lq — for 
even wave functions and Kq = and Lq = 1 for odd wave 
functions. If the Hamiltonian has no special symmetry, 
as for instance in the case of an asymmetric double- well, 
one could leave Lq indeterminate and fix it at the same 
time as E using conditions on the wave function or its 
derivative at two different points. These two conditions 
translate (in good approximation) into two polynomial 
equations in Lq and E and can be solved by Newton's 
method. 

The fact that Eqs. (8) determines entire functions pro- 
vided that V and G are polynomials can be inferred di- 
rectly from the fact that the coefficients will decrease as 
(n!)~ K for some positive power k to be determined and in 
general depending on the choice of G. As we will explain 
in more detail in Section IV, if the leading term in V is 
V21X 21 , one expects from Eq. (3) that for x large enough, 



<j>{x) ~ ±y / 2mV 2 ix l , 



and asymptotically, 



*(ar) oce-ornW 2 " 1 ^^ 1 



(9) 



(10) 



l+p—n 



Looking at the general expression for K given in Eq. (7), 
one sees that K will have the same asymptotic behav- 
ior provided that the integral of G grows not faster than 
x l+1 . If this is the case, then k = 1 / (Z + 1 ) . This behavior 
is well observed in empirical series. 

Note that if G grows faster than x , the coefficients 
decay more slowly and the procedure seem to be less ef- 
ficient. In the following, we will mostly discuss the case 
1 = 2. If we require that G is an odd polynomial growing 
not faster than x 2 , this means that G is homogeneous of 
degree 1. 



IV. QUANTIZATION FROM GLOBAL FLOW 
PROPERTIES 



In this section, we use the global properties of the flows 
associated with the Riccati equation to rephrase some im- 
plications of Sturm-Liouville theorem and to justify the 
asymptotic behavior given in Eq. (9). The main goal of 
this section is to provide a simple and intuitive picture of 
the bifurcation which occurs when the value of E is varied 
by a small amount above or below an energy eigenvalue. 
The main results of this section are summarized in Figs. 
1 and 2. 
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We consider the solutions of Eq. (3) obtained by vary- 
ing E with fixed initial values. It is convenient to intro- 
duce an additional parameter s and to rewrite the original 
equation as a 2-dimensional ODE with an s-independent 
r.h.s . 

K<j) = <j) 2 + 2m{E -V{x)) (11) 
x = l, (12) 

where the dot denotes the derivative with respect to s. 

The flows in the (x, (f) plane have some simple global 
properties that we now proceed to describe. We consider 
a solution (phase curve) with initial condition x = xq and 
4> = <f>o at s = 0. We assume that for these values the 
r.h.s of Eq. (11) is > 0. It will become clear later that if 
such a choice is impossible, a normalizable wave function 
cannot be constructed. With this assumption, the phase 
curve starts moving up and right as s increases, possibly 
going through simple poles with residues —h. This situ- 
ation persists unless the r.h.s. of (11) becomes zero. We 
call the separating curves defined by a zero for the r.h.s of 
Eq. (11), = ±y/2m(V(x)-E), ll WKB curves". Af- 
ter a phase curve crosses (horizontally) a WKB curve, 
it moves right and down. If it crosses the WKB curve 
again, we can repeat the discussion as at the beginning. 

At some point, we reach the "last" WKB curve (i.e., 
the farthest right). For x large enough, the potential is 
dominated by its largest power and the upper (lower) 
part of this last WKB curve has a strictly positive (nega- 
tive) slope. For such values of x, if a phase curve crosses 
the WKB curve, it will do so horizontally and move inside 
the region where the r.h.s. of Eq. (11) is negative. As s 
further increases, <f> decreases, but the phase curve can- 
not cross horizontally the lower part of the WKB curve 
which has a strictly negative slope. In the same region, if 
<j> has a pole, the curve reappears below the lower part of 
the WKB curve and will never take positive values again. 

In summary, if in the region described above, a phase 
curve crosses the WKB curve or develops a pole, then 
it cannot develop a pole again. The other logical pos- 
sibility is that the phase curve does none of the above. 
It is thus clear that for fixed E, we can always find a 
X such that if x > X, <j)(x) has no pole. Consequently 
the two terms involving <p m Eq. (3) cannot grow faster 
than 2m(E — V"). Otherwise, 2m(E — V) would become 
negligible and a pole would be necessary. At least one 
of these two terms needs to match 2m(E — V). Inspec- 
tion of the two possibilities leads to Eq. (9). Only the 
positive solution which follows asymptotically the upper 
WKB curve leads to a normalizable wave function. 

If we compare two phase curves with identical initial 
conditions but different E, the one with larger E initially 
lays above the other one. If the one with lower E has a 
first pole at x\, then the one with larger E has a first 
pole at some x < x\. Remembering that the poles of <p 
produce zeroes of "J, this rephrases the main idea behind 
the Sturm-Liouville theorem. An exact energy eigenstate 
E n is obtained when the wave function has its last zero 



at infinity. When E is fine-tuned to that value, <ft follows 
closely the upper branch of the WKB curve. This trajec- 
tory in unstable under small changes in E. If the energy 
is slightly increased with respect to E n , (f> develops a pole 
and reappears on the lower part of the WKB curve. If the 
energy is slightly decreased with respect to E n , (j> crosses 
the upper part of the WKB curve and reaches the lower 
part of the WKB curve. This is illustrated in Fig. 1 
in the case of the ground state of the quartic single-well 
anharmonic oscillator with m — 1/2, K = 1, Vi = 1 and 
Vi =0.1. All the figures in this section and the next two 
sections have been done with this particular example. 



H=x~2+p"2+0 . lx A 4 




1 2 3 4 5 6 



FIG. 1. Bifurcations of 4>(x) from the upper part of the 
WKB curve associated with the ground state energy Eo for 
energies E ± 1(T 5 , E Q ± 1CT 10 , E ± 1CT 15 , E Q ± 1(T 20 and 
E ± 1(T 25 (from left to right). 

The sensitive dependence on E is also present in the 
asymptotic behavior of K. If the energy is slightly in- 
creased with respect to E n , K reaches zero at a finite 
value of x. If the energy is slightly decreased with re- 
spect to E n , K increases rapidly. This is illustrated in 
Fig. 2 for the same example as in Fig. 1. 



H=x~2+p"2+0 . lx"4 




6 6.1 6.2 6.3 6.4 6.5 



FIG. 2. Bifurcations of K(x) from its trajectory for 
E = -B . The changes in E are ±1(T 30 , ±2 x 1(T 30 ,..., 
±1(T 29 

We now discuss the initial value (f> . For parity in- 
variant potentials, one only needs to consider the cases 
<t>o = (even or (j) = — oo (odd ^) at x = 0. For po- 
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tentials with no reflection symmetry, one needs to insure 
that the appropriate behavior is reached when x — > — oo. 
This can be implemented in good approximation by re- 
quiring that the wave function has also a zero at some 
large negative value x m i n . For potentials with a reflec- 
tion symmetry about another point x\ than the origin, 
one can impose that the wave function (K(x\) = 0) or 
its derivative (L(xi) — 0) vanish at that point. In all 
cases, we have an independent condition which allows to 
determine <fio- 

In summary, for x max large enough, the condition 

K(x max ) = (13) 

provides sharp upper bound on the energy levels. The 
lower part of Fig. 2 makes clear that as x max increases, 
sharper bounds are reached. For potentials that are 
not parity invariant, an additional condition has to be 
imposed. In all cases, one obtains polynomial equa- 
tions which can be solved for the energy levels given 
the potential or vice-versa using Newton's method. Note 
also that a sharp lower bound can be found by solving 
L{x max ) = 0. The fact that in Fig. 2, a zero of K at 
Eq + S corresponds to a zero of L at E — 6, suggests that 
the exact value should be very close to the average of the 
two bounds. 



V. G-DEPENDENCE OF THE BIFURCATION 

The strength of the bifurcation in K illustrated in Fig. 
2 can be approximately characterized by local exponents. 
If we consider the departure 5K(x) from K(x) calculated 
at some exact energy level E n , we expect the approximate 
behavior: 

SK(x) ~ C(E - E n )c xB . (14) 

In other words ln(|5iif (a;)|) is linear with a slope B inde- 
pendent of the choice of E and an intercept that varies 
like ln(|_E— E n \). This situation is approximately realized 
in the example considered before as shown in Fig. 3. We 
have checked in the same example that the sign of the 
energy difference plays no role. In other words, the same 
values of C and B can be used above and below E n . 
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6. 3 6.35 6.4 6. 45 6 .5 



FIG. 3. Natural logarithm of 8K(x) for E - E = 10 -30 
(lower set of point) and E — E — 1CF 28 (upper set of point). 
Lines are linear fits. 

The exponent B is not uniform. It increases with x 
and is G— dependent as shown in Fig. 4. The local val- 
ues of B have been calculated by fits in regions of width 
0.2 with central value displayed in the horizontal label of 
Fig. 4. 
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4 4.5 5 5.5 6 6.5 

XMAX 

FIG. 4. Value of B for various x and for G = —3x (empty 
hexagons), G = — 2x (filled squares) G = —x (crosses) and 
G = (triangles). The continuous lines have been drawn 
using Eq. (16). 

The change in B can be understood as follows. If E 
is changed from E n to E n + SE, then at some point we 
have a sudden transition from the upper to the lower 
WKB curve and asymptotically 

5*(x) oc SE c +TrTTW^ 2mV2 ' xl+1 . (15) 

Using Eq. (7) and expanding about x max , we obtain 
that, in good approximation, 



B ~ -(V2mV 2 ix l max - G{x max )) . (16) 

As shown in Fig. 4, this simple expression provides 
reasonable estimates of B. The slight underestimation 
comes in part from the fact that Eq. (16) does not take 
into account the harmonic term in V. Eq. (16) shows 
that we can increase the strength of the bifurcation near 
x m ax by increasing x max or -G(x max ). This allows us 
to "resolve" the energy more accurately. However, at 
the same time our numerical resolution of K(x max ) is af- 
fected and we need to take this effect into account. This 
question is treated in the next Section. In general, if we 
can establish that K(x max ) at an energy E very close 
to E n , can be calculated with some numerical accuracy 
6K num , we have the approximate numerical energy res- 
olution 

SE num - (x 5K num -e + ^(^^ 2mV2lX ™^ + $« ma:C dxG ^) 

(17) 
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VI. AN OPTIMAL CHOICE OF G 

In this Section, we show that from a numerical point of 
view, the choice of G is important. We discuss the ques- 
tion of an optimal choice, first with an example and then 
in general. We start with the calculation of the ground 
state in the case m — 1/2, h = 1,V% = 1 and V4 = 0.1. 
We discuss the estimation of the ground state energy us- 
ing the equation K(x max ) — with x max = 6. The fact 
that we use this finite value for x max creates an error in 
the 25-th digit (see Section VII). 

From the discussion of Section III, it is reasonable to 
limit the discussion to a gauge function of the form 

G{x) = -ax , (18) 

which using Eq. (7) implies that 

K(x) oc V(x)e^ ax2 . (19) 

With this restriction, the optimization problem is re- 
duced to the determination of a. As a increases through 
positive values, the features of '5 are exponentially am- 
plified, making the bifurcation displayed in Fig. 2 more 
violent. Ideally, we would like to take a as large as pos- 
sible. However, if a is too large, we may need too many 
coefficients K n to get a good approximation. If we con- 
sider the problem at a given order, the two requirements 
of sensitivity and accuracy result in a compromise which 
determines the optimal value of a. 

As explained in Section III, the choice of Eq. (18) 
guarantees a suppression of the form (nl)~ for the coef- 
ficients of L and K. However, the choice of a still affects 
significantly the behavior of these coefficients as shown 
in Fig. 5. 
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25 50 75 100 125 150 175 200 

n 

FIG. 5. ln(|is: n 6 n |) versus n, for G = (triangles), G = -x 
(filled squares), G = — 2x (crosses) and G = —3x (empty 
squares) . 

The quantity KnX^ax i s relevant to decide at which 
order we need to truncate the series in order to get a 
good estimate of K{x max ). For instance, if we require 
knowing K{x max ) with errors of order 1, we need about 
100 coefficients for a = 2 but more than 150 for a = 0. 
The corresponding values for a = 1 and 3 fall between 



these two values, indicating that a = 2 is close to optimal. 
This estimate is confirmed by an analysis of the depen- 
dence of K n on a. Sample values are shown in Fig. 6. 
We observe rapid oscillations (that we will not attempt 
to explain) and slowly varying amplitudes which have a 
minimum slightly below 2. Note that on the logarithmic 
scale of Fig. 6, the zeroes of K n give —00, however due 
to the discrete sampling of a, it just generates isolated 
dots on the graphs. Note also that in Figs. 5 and 6, the 
coefficients have been calculated for an an accurate value 
of the ground state energy. 
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FIG. 6. ln(|if n |) versus a, for n=60 (upper set), 70 (next 
set), 80 (next set) and 90 (lower set). 

The behavior of the K n calculated at value of E suffi- 
ciently close to an eigenvalue, can be understood by using 
the asymptotic form 

K(x max ) oc e ^-Th^^L-r- . (20) 

We emphasize that the relative sign between the two 
terms in the exponential is opposite than in Eq. (17), 
because we are now on the upper WKB curve. For a = 0, 
Eq. (20) provides a rough estimate of KnX 7 ^^. Remem- 
bering the minus sign in the parametrization of G (Eq. 
(18)), we see that if a is given a small positive value, the 
argument of the exponential in Eq. (20) decreases and 
we can obtain comparable accuracy with less terms in 
the expansion. Naively, our optimum choice is obtained 
when the two terms in the exponential cancel. In the 
general case, this amounts to having 

^2mV 2l x l +l x ~ -(/ + 1) / dxG(x) . (21) 

Jo 

For the particular example considered here, this cancel- 
lation is obtained for a = (2/3)V0Ax max ~ 1.27. It is 
clear that when the two terms cancel, subleading terms 
neglected in Eq. (9) should be taken into account. How- 
ever, in several examples, we found that this simple pro- 
cedure gives results close to what is found empirically. 

We now address the more general question of determin- 
ing the G-dependence of the number of significant digits 
that can be obtained from the condition K{x max ) = 
using an expansion of K truncated at a given order. For 
the example considered before in this section, we see from 



G 



Fig. 7 that, for instance for a truncation at order 100, the 
most accurate answer is obtained for a ~ 1.6. It is worth 
noting that for this value of a, one gains more than 15 
significant digits compared to the G = case! This figure 
also indicates, that as expected, the best possible answer 
(in the present case, 25 significant digits) can always be 
achieved by calculating enough coefficients. 
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FIG. 7. Number of significant digits for Eq versus a using 
the condition K(6) = with expansions of order 50 (empty 
diamonds), 75 (filled squares), 100 (crosses), 125 (empty 
squares) and 150 (stars). 

Using Eq. (17) and Fig. (6), we were able to reproduce 
approximately the left part of Fig. 7 (0 < a < 1). To 
give a specific example, at order 100, when one changes 
a from to 1, SK num becomes 4 orders of magnitude 
smaller and the factor e"^) 1 ™.. improves the resolution 
by almost 8 orders of magnitude. This approximately ac- 
counts for the gain of 11 significant digits observed in Fig. 
7. A detailed understanding of the figure in the region 
1 < a < 2 is beyond what can be accomplished using 
the asymptotic form of the wave function. However, the 
naive estimate of Eq. (21) provides a reasonable estimate 
of the location of the optimal a. 

It should be noted that an ansatz of the form of Eq. 
(19) with a = 1 has been used in Rcf. [12] and that 
the fact that varying a could improve the numerical effi- 
ciency was found empirically in Ref. [13]. Eq. (21) can 
be used to understand these results. For instance, for 
H = p 2 + x 2 + x 8 , we can obtain a very accurate result 
with x max = 2.8 (see Section VII). According to Eq. (21) 
the optimal value of a in this case is a = (2/5)x 3 nax ~ 8.8 
which is slightly below the value (m 10) suggested in [13]. 
Note also that equivalently good results can be obtained 
using G = —bx 3 . 



VII. APPROXIMATE ERROR FORMULAS 

In this Section, we discuss the intrinsic error SE = 
E(x max ) - E(oo) where E(x max ) is defined by 
ip(x max ,E(x max )) = 0, for a given energy level. We 
emphasize that SE is the error due to the finiteness of 
Xmax independently of practical considerations regard- 



ing the numerical estimation of E(x max ) which is as- 
sumed to be known with an error much smaller than 
5E in this section. We use the familiar parametrization 
of the quadratic term of the potential, Vi — \muj 2 and 
we restore the dependence on Ti and m. The error for 
the ground state of the harmonic oscillator has been esti- 
mated in Eq. (4) of Ref. [20] . Using the asymptotic form 
of the integral in this equation, we obtain 

SE harm. _ g (_0) ± ^ So /h _ ^) 

with 

f + °° 1 

S = dt-m{{x c {t)f + u (x c (t)) ) = mujx 2 max 

(23) 

and x c {t) = ^maxe - ' 1 ''* - * ' . This corresponds to semi- 
classical approximation where the contribution of the 
large field configurations are obtained by calculating the 
quadratic fluctuations with respect to x c (t). The anhar- 
monic corrections can be approximated to lowest order 
in the the anharmonic couplings by adding a term S an h 
to So in the exponent of Eq. (22) with 

/+00 
dtV anh (x c (t)) , (24) 
-00 

and Vanh is the anharmonic part of the potential. Our 
final pcrturbative estimate is thus 

5E ~ 5££«™. e -£ 3 . =2 (jk)v-2^ ox ( 25 j 

This estimate is accurate if the V^j and x max are small 
enough. We expect that for the excited states, approx- 
imate formulas of the form of Eq. (25) multiplied by a 
polynomial should hold. 

When A or x max become too large, Eq. (25) is not 
adequate. To obtain a better approximation, we use 

d 

a i>(x max ,E(x m ax)) = , (26) 

and the asymptotic behavior of We estimate that 
d^> I dE is of the order of the non-normalizable WKB so- 
lution and as a consequence, SE has the asymptotic form 

SE ~ P(x 

in a. r mXmax)) 2 , (27) 

where P is a polynomial. This form is correct for the 
ground state of the harmonic oscillator. In the case where 
the leading term of V is VnX 1 , this implies the asymp- 
totic order of magnitude estimate 

SE w e^W^ 1 ™ , (28) 

We have tested the two approximate errors formulas 
given above (Eqs. (25) and (28)) for the ground state 
corresponding to V an h = Ax 4 . We used the numerical 
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values h = m = oj = 1 and A = 0.1 The results are 
shown in Fig. 8. We see that for small values of x maxi 
the perturbative estimate of Eq. (25) corrects properly 
the harmonic result. However when the 
Eq. (28) gives better results. If the left part of the graph 
is displayed with a log-log scale, it is approximately linear 
with a slope close to 3. In Fig. (8), the proportionality 
constant not given by Eq. (28) has been determined by 
fitting the 5 last data points on the right of the figure. 
We conclude that by combining the two approximations 
it is possible to get a reasonable estimate of the errors on 
E over a wide range of x max . 



LAMBDA =0.1 




3 4 5 6 7 

XMAX 

FIG. 8. ln( SEo) as a function of Xmax for A=0.1 (black 
dots). The continuous lines are from top to bottom on the 
left of the figure: the harmonic case (Eq. (22)), Eq. (25) with 
Va = 0.1 (fits the dots well on the left of the figure), Eq. (28) 
(fits the dots well on the right of the figure). 

We have tested Eq. (28) for other potentials. For in- 
stance, for H = p 2 + x 2 + x 8 , in order to get 30 significant 
digit, we estimated that x max ~ 2.8. We found that the 
difference between the ground state energy found from 
the condition K = (upper bound) and L = (lower 
bound) differed in the 30th significant digits. 



VIII. A CHALLENGING TEST 

The only practical limitation of the method proposed 
here is that in some cases the relevant details of the po- 
tential appear in widely separated regions, forcing us to 
calculate a huge number of coefficients with many signif- 
icant digits. A simple example where such problem may 
occur is the symmetric double-well with a small quartic 
coupling where the separation between the wells goes like 
the inverse square root of the quartic coupling. 

In Ref. [23], the lowest even and odd energies were 
calculated for a potential with m = 1, K = 1, Vi = 
-1/4, V 4 = 1/2000 with 180 significant digits. Remark- 
ably, the authors were able to reproduce the 110 sig- 
nificant digits of the splitting between these two states 
by calculating instanton effects. We have reproduced the 
180 digits of both states using an expansion of order 1700 
for K and a value of x max = 46. The calculations were 



performed with 700 digit arithmetic. The calculation of 
one level with such a procedure takes less than two hours 
with MATHEMATICA on an unexpensive laptop using 
Pentium3. The computation time increases with the ac- 
curacy required. In order to fix the ideas, it takes less 
than 2 minutes minutes to reproduce the first 120 digits 
in the above calculation. 



IX. THE MULTIVARIABLE CASE 

The basic equations presented in Section II can be ex- 
tended when the single variable x is replaced by a N- 
dimensional vector x. In Eq. (2), becomes a vector 
<fr and the integral a line integral. In order to guarantee 
that the wave function is independent of the choice of 
the line, we require that the curl of <p vanishes. Eq. (3) 
becomes: 

hV$=$-$+2m(E-V) . (29) 

Using <f> = L/K, we write as previously 

hVL + 2m(V - E)K + G ■ L = (30) 
TiVK + L + GK = , (31) 

with G(x) unspecified at this point. These equations im- 
ply the multivariable Riccati equation (29) multiplied by 
K 2 . Near a zero of K, these equations imply the same 
singularity as Eq. (29). After using Eq. (31), the condi- 
tion that (f> has no curl reads 

VjLy) + GtyL^ = VjLi + G {j )L {i) . (32) 

The parenthesis for the vector indices are used in order to 
distinguish these indices from the order in a power series 
expansion used later. 

The transformation Eqs. (7) can vectorized trivially 
with Q treated as a scalar. In the expression of K given 
by Eq. (7), the integral becomes a line integral and we 
require that G(x) has a vanishing curl. This condition 
is also necessary to establish that different derivatives 
acting on K commute. 

The choice of coordinates to be used depends on the 
choice of boundary conditions imposed. If we require 

to vanish on a large hypersphere, hypcrspherical har- 
monics should be used. If we require \I/ to vanish on hy- 
percubes (as suggested for lattice problems in Ref. [20]) 
cartesian coordinates should be used. To fix the ideas, 
let us consider the case of cartesian coordinates for two 
variables Xi and X2 with boundary conditions on a rect- 
angle. We expand K(x 1 ,x 2 ) = J2 m , n >o K rn, n xT l x'2 and 

similar expansions for the two components of L. The 
coefficients can be constructed order by order, with the 
order of K m , n defined as m + n. The terms with one 
derivative yield the higher order terms. For instance, for 
K, we obtain equations providing K m +i^ n and K m n ^-i 
in terms of coefficients of higher order just as in Eq. (8). 



8 



A detailed construction shows that if V(x\,x 2 ) has no 
special symmetry, we can determine all the coefficient 
up to a given order / provided that we supply the val- 
ues of two coefficients at each intermediate order (for in- 
stance L m .o for m < I). These coefficients together with 
E are fixed by the boundary conditions K(x\ m i n , x 2 ) — 
K(x lmax ,x 2 ) = K(xi,x 2m in) = K(x 1 ,x 2ma x) = 0. Tak- 
ing derivatives with respect to the free variables x\ and 
x 2 , and setting these variables to 0, we obtain an infinite 
set of conditions. The truncation of this set, together 
with the truncation of the expansion in the other vari- 
able must be studied carefully. If we consider the special 
case where the problem can be solved by separation of 
variables, we see that it is important to maintain a uni- 
form accuracy for all the conditions. If all the coefficients 
have been calculated up to order /, this can achieved in 
the following way. We retain of the order of 1/2 deriva- 
tives of the four conditions in such a way that we get 
exactly 21 + 3 conditions which can be expanded up to 
an order close to 1/2 in the remaining variable. A prac- 
tical implementation of this program is in progress. 



ables. It remains to be determined if the simultaneous 
solution of many polynomial equations can be accom- 
plished with a reasonable accuracy. For these problems, 
the fact that a judicious choice of the arbitrary functions 
G allows to decrease the order of the expansions may be 
crucial. 
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X. CONCLUSIONS 



In conclusion, we have shown that accurate estimates 
of the energy levels of arbitrary polynomial potentials 
bounded from below can be obtained by solving poly- 
nomial equations. The fact that the function L and K 
are entire guarantees that if we calculate enough terms 
we will gain proper control of the asymptotic behavior of 
the wave function. Reaching this goal is in general a dif- 
ficult task which often requires guesswork and analytical 
continuations (see e.g., Ref. [24]). Here, the convergence 
of the procedure is guaranteed and the order at which 
we can terminate the expansion in order to reach a given 
accuracy can be estimated. In addition, a systematic un- 
derstanding and control of the errors due to the finite 
value of x max has been achieved. 

The understanding of the gauge invariance of the ba- 
sic equations proposed here completely resolves the issues 
raised from our initial proposal [15]. By varying G, from 
to -(f), we can interpolate between a situation where K 
is the wave function to another situation where K = 1 
and L = <j>. However, for every other choice of G, only 
the ratio L/K has a direct physical meaning. By prop- 
erly chosing G, we can at the same time improve the 
convergence of K and amplify the bifurcation toward the 
the non-normalizable behavior. 

The extreme accuracy obtained for two widely sepa- 
rated wells indicates that for reasonably complicated po- 
tential, the number of terms that needs to be calculated is 
not prohibitive. We intend to use this method to test an- 
alytical results regarding the role of large configurations 
in the path-integral and to test semi-classical treatment 
of potentials with asymmetric wells [6,7]. 

The method can be extended in the case of several vari- 
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